main_pipeline
wentao
2023/01/28
写在前面
2023年01月28日，更新ggClusterNet文档。
R包安装
导入R包
#--导入所需R包#-------
library(phyloseq)
library(igraph)
library(network)
library(sna)
library(tidyverse)
library(ggClusterNet)输入方式
输入方式一
直接输入phyloseq格式的数据。
#-----导入数据#-------
data(ps)输入方式二
可以从https://github.com/taowenmicro/R-_function下载数据，构造phylsoeq文件。自己的数据也按照网址示例数据进行准备。虽然phylsoeq对象不易用常规手段处理，但是组学数据由于数据量比较大，数据注释内容多样化，所以往往使用诸如phyloseq这类对象进行处理，并简化数据处理过程。ggClusterNet同样使用了phyloseq对象作为微生物网络的分析。
phyloseq对象构建过程如下，网络分析主要用到otu表格，后续pipeline流程可能用到分组文件metadata，如果按照分类水平山色或者区分模块则需要taxonomy。这几个部分并不是都必须加入phyloseq对象中，可以用到那个加那个。
或者直接从网站读取，但是由于github在国外，所以容易失败
使用ggCLusterNet进行网络分析的过程
corMicro函数用于计算网络相关
按照丰度过滤微生物表格，并却计算相关矩阵，按照指定的阈值挑选矩阵中展示的数值。调用了psych包中的corr.test函数，使用三种相关方法。 N参数提取丰度最高的150个OTU； method.scale参数确定微生物组数据的标准化方式，这里我们选用TMM方法标准化微生物数据。
#-提取丰度最高的指定数量的otu进行构建网络
#----------计算相关#----
result = corMicro(ps = ps,
N = 150,
method.scale = "TMM",
r.threshold=0.8,
p.threshold=0.05,
method = "spearman"
)
#--提取相关矩阵
cor = result[[1]]
cor[1:6,1:6]## ASV_1 ASV_28 ASV_125 ASV_135 ASV_8 ASV_106
## ASV_1 1 0 0 0 0 0
## ASV_28 0 1 0 0 0 0
## ASV_125 0 0 1 0 0 0
## ASV_135 0 0 0 1 0 0
## ASV_8 0 0 0 0 1 0
## ASV_106 0 0 0 0 0 1#-网络中包含的OTU的phyloseq文件提取
ps_net = result[[3]]
#-导出otu表格
otu_table = ps_net %>%
vegan_otu() %>%
t() %>%
as.data.frame()制作分组,我们模拟五个分组
这是网络布局的基础，无论是什么聚类布局，都需要制作一个分组文件，这个文件有两列，一列是节点，一列是分组信息，这个分组信息名称为：group。这个文件信息就是用于对节点进行分组，然后按照分组对节点归类，使用包中可视化函数计算节点位置。
注意分组文件的格式，分为两列，第一列是网络中包含的OTU的名字，第二列是分组信息，同样的分组标记同样的字符。
#--人工构造分组信息：将网络中全部OTU分为五个部分，等分
netClu = data.frame(ID = row.names(otu_table),group =rep(1:5,length(row.names(otu_table)))[1:length(row.names(otu_table))] )
netClu$group = as.factor(netClu$group)
head(netClu)## ID group
## 1 ASV_1 1
## 2 ASV_28 2
## 3 ASV_125 3
## 4 ASV_135 4
## 5 ASV_8 5
## 6 ASV_106 1PolygonClusterG 根据分组，计算布局位置坐标
不同的模块按照分组聚集成不同的圆，并且圆形的大小一样。如果一个分组只有一个点，则这个点坐落在圆心位置。
#--------计算布局#---------
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
head(node)## X1 X2 elements
## ASV_1 0.000000e+00 10.5 ASV_1
## ASV_14 2.598076e+00 9.0 ASV_14
## ASV_72 2.598076e+00 6.0 ASV_72
## ASV_109 3.673819e-16 4.5 ASV_109
## ASV_71 -2.598076e+00 6.0 ASV_71
## ASV_88 -2.598076e+00 9.0 ASV_88nodeadd 节点注释的：用otu表格和分组文件进行注释
nodeadd函数只是提供了简单的用注释函数，用户可以自己在node的表格后面添加各种注释信息。
tax_table = ps_net %>%
vegan_tax() %>%
as.data.frame()
#---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
nodes[1:6,1:6]## X1 X2 elements Kingdom Phylum
## ASV_1 0.000000 10.5000000 ASV_1 Bacteria Actinobacteria
## ASV_10 3.792382 -8.9964485 ASV_10 Bacteria Proteobacteria
## ASV_100 7.269286 -5.1349547 ASV_100 Bacteria Proteobacteria
## ASV_101 5.911236 5.0628075 ASV_101 Bacteria Proteobacteria
## ASV_102 4.278276 1.3951201 ASV_102 Bacteria Proteobacteria
## ASV_103 -8.359017 -0.4411929 ASV_103 Bacteria Proteobacteria
## Class
## ASV_1 Actinobacteria
## ASV_10 Alphaproteobacteria
## ASV_100 Betaproteobacteria
## ASV_101 Unassigned
## ASV_102 Gammaproteobacteria
## ASV_103 Alphaproteobacteria计算边
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
head(edge)## X2 Y2 OTU_2 OTU_1 weight X1 Y1 cor
## 1 2.598076 9.0000000 ASV_14 ASV_172 -0.8020673 4.278276 3.2492221 -
## 2 7.131446 5.3221711 ASV_28 ASV_20 0.8041345 7.755181 -0.6122717 +
## 3 4.533370 0.8221711 ASV_140 ASV_64 0.8484250 1.818041 -4.5620057 +
## 4 7.014193 -4.5620057 ASV_166 ASV_43 -0.8064018 -2.853170 8.4270510 -
## 5 7.014193 -4.5620057 ASV_166 ASV_172 -0.8016529 4.278276 3.2492221 -
## 6 1.818041 -7.5620057 ASV_29 ASV_43 -0.8415076 -2.853170 8.4270510 -出图
p <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
theme(panel.background = element_blank()) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p
# ggsave("cs1.pdf",p,width = 8,height = 6)模拟不同的分组–可视化
模拟不同分组效果展示：3个分组
这是网络布局的基础，无论是什么聚类布局，都需要制作一个分组文件，这个文件有两列，一列是节点，一列是分组信息，这个分组信息名称必须设定为：group。
netClu = data.frame(ID = row.names(tax_table),group =rep(1,length(row.names(tax_table)))[1:length(row.names(tax_table))] )
netClu$group = as.factor(netClu$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs2.pdf",pnet,width = 7,height = 5.5)模拟不同的分组查看效果：8个分组
netClu = data.frame(ID = row.names(cor),group =rep(1:8,length(row.names(cor)))[1:length(row.names(cor))] )
netClu$group = as.factor(netClu$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs3.pdf",pnet,width = 6,height = 5)微生物分类分组可视化
#----------计算相关#----
result = corMicro (ps = ps,
N = 200,
method.scale = "TMM",
r.threshold=0.8,
p.threshold=0.05,
method = "spearman"
)
#--提取相关矩阵
cor = result[[1]]
# head(cor)
#-网络中包含的OTU的phyloseq文件提取
ps_net = result[[3]]
#-导出otu表格
otu_table = ps_net %>%
vegan_otu() %>%
t() %>%
as.data.frame()
tax = ps_net %>% vegan_tax() %>%
as.data.frame()
tax$filed = tax$Phylum
group2 <- data.frame(ID = row.names(tax),group = tax$Phylum)
group2$group =as.factor(group2$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs3.pdf",pnet,width = 6,height = 5)微生物分类可视化布局优化1-圆环半径调整PolygonRrClusterG
结果发现这些高丰度OTU大部分属于放线菌门和变形菌门，其他比较少。所以下面我们按照OTU数量的多少，对每个模块的大小进行重新调整。
result2 = PolygonRrClusterG(cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs4.pdf",pnet,width = 6,height = 5)微生物分类可视化布局优化2-model_filled_circle
用实心点作为每个模块的布局方式
set.seed(12)
#-实心圆2
result2 = model_filled_circle(cor = cor,
culxy =TRUE,
da = NULL,# 数据框，包含x,和y列
nodeGroup = group2,
mi.size = 1,# 最小圆圈的半径，越大半径越大
zoom = 0.3# 不同模块之间距离
)## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs5.pdf",pnet,width = 6,height = 5)微生物分类可视化布局优化3 model_maptree_group
用实心点作为每个模块布局方式2：model_maptree_group，智能布局不同分组之间的距离，在美学上特征更明显一点。
set.seed(12)
#-实心圆2
result2 = model_maptree_group(cor = cor,
nodeGroup = group2,
)
# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs6.pdf",pnet,width = 6,height = 5)按照网络模块化分组
模块布局算法 model_maptree_group
按照物种组成分类完成网络分析其实并不常用，更多的是按照模块分组，进行网络可视化。
#--modulGroup函数用于计算模块并整理成分组信息
netClu = modulGroup( cor = cor,cut = NULL,method = "cluster_fast_greedy" )
result2 = model_maptree_group(cor = cor,
nodeGroup = group2,
)
# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
# head(nodes)
nodes2 = nodes %>% inner_join(netClu,by = c("elements" = "ID"))
nodes2$group = paste("Model_",nodes2$group,sep = "")
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = group,size = mean),pch = 21, data = nodes2) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs7.pdf",pnet,width = 6,height = 5)模块布局算法 model_maptree2
使用升级的model_maptree2：不在可以将每个模块独立区分，而是将模块聚拢，并在整体布局上将离散的点同这些模块一同绘制到同心圆内。控制了图形的整体布局为圆形。
result2 = model_maptree2(cor = cor,
method = "cluster_fast_greedy"
)
# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)## X1 X2 elements Kingdom Phylum
## ASV_1 -13.8901865 -0.03852081 ASV_1 Bacteria Actinobacteria
## ASV_10 -0.4097123 3.33964715 ASV_10 Bacteria Proteobacteria
## ASV_100 -1.1227948 -12.61245151 ASV_100 Bacteria Proteobacteria
## ASV_101 5.7161404 20.40310410 ASV_101 Bacteria Proteobacteria
## ASV_102 -18.4899174 -10.64072413 ASV_102 Bacteria Proteobacteria
## ASV_103 16.2557654 10.33992920 ASV_103 Bacteria Proteobacteria
## Class Order Family Genus
## ASV_1 Actinobacteria Actinomycetales Thermomonosporaceae Unassigned
## ASV_10 Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium
## ASV_100 Betaproteobacteria Burkholderiales Comamonadaceae Hydrogenophaga
## ASV_101 Unassigned Unassigned Unassigned Unassigned
## ASV_102 Gammaproteobacteria Unassigned Unassigned Unassigned
## ASV_103 Alphaproteobacteria Rhizobiales Phyllobacteriaceae Phyllobacterium
## Species mean
## ASV_1 Unassigned 1473.61111
## ASV_10 Unassigned 428.27778
## ASV_100 Hydrogenophaga_intermedia 61.00000
## ASV_101 Unassigned 60.72222
## ASV_102 Unassigned 64.00000
## ASV_103 Phyllobacterium_bourgognense 53.77778nodes2 = nodes %>% inner_join(netClu,by = c("elements" = "ID"))
nodes2$group = paste("Model_",nodes2$group,sep = "")
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
data = edge, size = 0.5) +
geom_point(aes(X1, X2,fill = group,size = mean),pch = 21, data = nodes2) +
scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# labs( title = paste(layout,"network",sep = "_"))+
# geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
# theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
# ggsave("cs8.pdf",pnet,width = 6,height = 5)上千OTU相关性计算测试model_igraph2布局
这个布局最近几年文章上使用非常多。
# cor_Big_micro2 增加了标准化方法和p值矫正方法
result = cor_Big_micro2(ps = ps,
N = 1000,
r.threshold=0.85,
p.threshold=0.05,
method = "pearson",
scale = FALSE
)
#--提取相关矩阵
cor = result[[1]]
dim(cor)## [1] 1000 1000# model_igraph2
result2 <- model_igraph2(cor = cor,
method = "cluster_fast_greedy",
seed = 12
)
node = result2[[1]]
dim(node)## [1] 293 3dat = result2[[2]]
head(dat)## orig_model model color OTU X1 X2
## 1 45 mini_model #C1C1C1 ASV_1591 7.835825 12.339063
## 2 15 model_15 #50A9AF ASV_688 -1.523618 -9.873582
## 3 34 model_34 #FEF6B1 ASV_8 4.877237 13.710152
## 4 6 model_6 #F99655 ASV_174 2.811909 -10.288624
## 5 6 model_6 #F99655 ASV_716 4.314573 -10.310663
## 6 34 model_34 #FEF6B1 ASV_106 3.768050 13.839837tem = data.frame(mod = dat$model,col = dat$color) %>%
dplyr::distinct( mod, .keep_all = TRUE)
col = tem$col
names(col) = tem$mod
#---node节点注释#-----------
otu_table = as.data.frame(t(vegan_otu(ps)))
tax_table = as.data.frame(vegan_tax(ps))
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)## X1 X2 elements Kingdom Phylum
## ASV_10 -2.259642 10.577060 ASV_10 Bacteria Proteobacteria
## ASV_100 -7.994668 8.125110 ASV_100 Bacteria Proteobacteria
## ASV_102 6.730077 -13.330438 ASV_102 Bacteria Proteobacteria
## ASV_104 -2.682913 12.612126 ASV_104 Bacteria Proteobacteria
## ASV_1043 3.463727 -9.298511 ASV_1043 Bacteria Acidobacteria
## ASV_1048 12.307434 -7.316681 ASV_1048 Bacteria Actinobacteria
## Class Order Family Genus
## ASV_10 Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium
## ASV_100 Betaproteobacteria Burkholderiales Comamonadaceae Hydrogenophaga
## ASV_102 Gammaproteobacteria Unassigned Unassigned Unassigned
## ASV_104 Betaproteobacteria Burkholderiales Comamonadaceae Polaromonas
## ASV_1043 Acidobacteria_Gp10 Unassigned Unassigned Gp10
## ASV_1048 Actinobacteria Actinomycetales Microbacteriaceae Agromyces
## Species mean
## ASV_10 Unassigned 428.277778
## ASV_100 Hydrogenophaga_intermedia 61.000000
## ASV_102 Unassigned 64.000000
## ASV_104 Unassigned 60.444444
## ASV_1043 Unassigned 4.833333
## ASV_1048 Agromyces_indicus 4.722222#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
colnames(edge)[8] = "cor"
head(edge)## X2 Y2 OTU_2 OTU_1 weight X1 Y1 cor
## 1 7.835825 12.339063 ASV_1591 ASV_710 0.9060991 7.958023 13.239291 +
## 2 -1.523618 -9.873582 ASV_688 ASV_596 0.9262611 -2.153803 -10.584342 +
## 3 -1.523618 -9.873582 ASV_688 ASV_780 0.8884699 -1.365654 -10.865238 +
## 4 -1.523618 -9.873582 ASV_688 ASV_20 0.8557562 -1.207825 -8.872789 +
## 5 4.877237 13.710152 ASV_8 ASV_106 0.8734179 3.768050 13.839837 +
## 6 4.877237 13.710152 ASV_8 ASV_334 0.9389280 4.640294 14.523675 +tem2 = dat %>%
dplyr::select(OTU,model,color) %>%
dplyr::right_join(edge,by =c("OTU" = "OTU_1" ) ) %>%
dplyr::rename(OTU_1 = OTU,model1 = model,color1 = color)
head(tem2)## OTU_1 model1 color1 X2 Y2 OTU_2 weight X1
## 1 ASV_1591 mini_model #C1C1C1 7.958023 13.239291 ASV_710 0.9060991 7.835825
## 2 ASV_688 model_15 #50A9AF -2.153803 -10.584342 ASV_596 0.9262611 -1.523618
## 3 ASV_688 model_15 #50A9AF -1.365654 -10.865238 ASV_780 0.8884699 -1.523618
## 4 ASV_688 model_15 #50A9AF -1.207825 -8.872789 ASV_20 0.8557562 -1.523618
## 5 ASV_8 model_34 #FEF6B1 3.768050 13.839837 ASV_106 0.8734179 4.877237
## 6 ASV_8 model_34 #FEF6B1 4.640294 14.523675 ASV_334 0.9389280 4.877237
## Y1 cor
## 1 12.339063 +
## 2 -9.873582 +
## 3 -9.873582 +
## 4 -9.873582 +
## 5 13.710152 +
## 6 13.710152 +tem3 = dat %>%
dplyr::select(OTU,model,color) %>%
dplyr::right_join(edge,by =c("OTU" = "OTU_2" ) ) %>%
dplyr::rename(OTU_2 = OTU,model2 = model,color2 = color)
head(tem3)## OTU_2 model2 color2 X2 Y2 OTU_1 weight X1
## 1 ASV_1591 mini_model #C1C1C1 7.835825 12.339063 ASV_710 0.9060991 7.958023
## 2 ASV_688 model_15 #50A9AF -1.523618 -9.873582 ASV_596 0.9262611 -2.153803
## 3 ASV_688 model_15 #50A9AF -1.523618 -9.873582 ASV_780 0.8884699 -1.365654
## 4 ASV_688 model_15 #50A9AF -1.523618 -9.873582 ASV_20 0.8557562 -1.207825
## 5 ASV_8 model_34 #FEF6B1 4.877237 13.710152 ASV_106 0.8734179 3.768050
## 6 ASV_8 model_34 #FEF6B1 4.877237 13.710152 ASV_334 0.9389280 4.640294
## Y1 cor
## 1 13.239291 +
## 2 -10.584342 +
## 3 -10.865238 +
## 4 -8.872789 +
## 5 13.839837 +
## 6 14.523675 +tem4 = tem2 %>%inner_join(tem3)
head(tem4)## OTU_1 model1 color1 X2 Y2 OTU_2 weight X1
## 1 ASV_1591 mini_model #C1C1C1 7.958023 13.239291 ASV_710 0.9060991 7.835825
## 2 ASV_688 model_15 #50A9AF -2.153803 -10.584342 ASV_596 0.9262611 -1.523618
## 3 ASV_688 model_15 #50A9AF -1.365654 -10.865238 ASV_780 0.8884699 -1.523618
## 4 ASV_688 model_15 #50A9AF -1.207825 -8.872789 ASV_20 0.8557562 -1.523618
## 5 ASV_8 model_34 #FEF6B1 3.768050 13.839837 ASV_106 0.8734179 4.877237
## 6 ASV_8 model_34 #FEF6B1 4.640294 14.523675 ASV_334 0.9389280 4.877237
## Y1 cor model2 color2
## 1 12.339063 + mini_model #C1C1C1
## 2 -9.873582 + model_15 #50A9AF
## 3 -9.873582 + model_15 #50A9AF
## 4 -9.873582 + model_15 #50A9AF
## 5 13.710152 + model_34 #FEF6B1
## 6 13.710152 + model_34 #FEF6B1edge2 = tem4 %>% mutate(color = ifelse(model1 == model2,as.character(model1),"across"),
manual = ifelse(model1 == model2,as.character(color1),"#C1C1C1")
)
col_edge = edge2 %>% dplyr::distinct(color, .keep_all = TRUE) %>%
select(color,manual)
col0 = col_edge$manual
names(col0) = col_edge$color
library(ggnewscale)
p1 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = color),
data = edge2, size = 1) +
scale_colour_manual(values = col0)
# ggsave("./cs1.pdf",p1,width = 16,height = 14)
p2 = p1 +
new_scale_color() +
geom_point(aes(X1, X2,color =model), data = dat,size = 4) +
scale_colour_manual(values = col) +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
theme(panel.background = element_blank()) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p2
# ggsave("./cs2.pdf",p2,width = 12,height = 10)网络属性和节点属性
网络性质计算
22年6月升级后版本包括了16项网络属性，包括周集中老师21年NCC文章中全部属性
dat = net_properties(igraph)
head(dat)## value
## num.edges 122.0000000
## num.pos.edges 101.0000000
## num.neg.edges 21.0000000
## num.vertices 49.0000000
## connectance 0.1037415
## average.degree 4.9795918# 升级后包含的网络属性更多
dat = net_properties.2(igraph,n.hub = T)
head(dat,n = 16)## value
## num.edges(L) 122.0000000
## num.pos.edges 101.0000000
## num.neg.edges 21.0000000
## num.vertices(n) 49.0000000
## Connectance(edge_density) 0.1037415
## average.degree(Average K) 4.9795918
## average.path.length 2.1710403
## diameter 4.7234210
## edge.connectivity 0.0000000
## mean.clustering.coefficient(Average.CC) 0.5690439
## no.clusters 5.0000000
## centralization.degree 0.1879252
## centralization.betweenness 0.1533413
## centralization.closeness 1.1468429
## RM(relative.modularity) 0.2573947
## the.number.of.keystone.nodes 1.0000000dat = net_properties.3(igraph,n.hub = T)
head(dat,n = 16)## value
## num.edges(L) 122.0000000
## num.pos.edges 101.0000000
## num.neg.edges 21.0000000
## num.vertices(n) 49.0000000
## Connectance(edge_density) 0.1037415
## average.degree(Average K) 4.9795918
## average.path.length 2.1710403
## diameter 4.7234210
## edge.connectivity 0.0000000
## mean.clustering.coefficient(Average.CC) 0.5690439
## no.clusters 5.0000000
## centralization.degree 0.1879252
## centralization.betweenness 0.1533413
## centralization.closeness 1.1468429
## RM(relative.modularity) -8.4867310
## the.number.of.keystone.nodes 1.0000000# 增加了网络模块性（modularity.net ）和随机网络模块性（modularity_random ）
# dat = net_properties.4(igraph,n.hub = F)
# head(dat,n = 16)节点性质计算
nodepro = node_properties(igraph)
head(nodepro)## igraph.degree igraph.closeness igraph.betweenness igraph.cen.degree
## ASV_28 3 0.2569444 2.333333 3
## ASV_48 1 0.2032967 0.000000 1
## ASV_34 5 0.3274336 140.216667 5
## ASV_18 4 0.2846154 20.916667 4
## ASV_20 3 0.2569444 36.000000 3
## ASV_44 11 0.3627451 7.874009 11Zipi基于模块对OTU进行分类
result = cor_Big_micro2(ps = ps,
N = 500,
r.threshold=0.6,
p.threshold=0.05,
# method = "pearson",
scale = FALSE
)
#--提取相关矩阵
cor = result[[1]]
result4 = nodeEdge(cor = cor)
#提取变文件
edge = result4[[1]]
#--提取节点文件
node = result4[[2]]
igraph = igraph::graph_from_data_frame(edge, directed = FALSE, vertices = node)
res = ZiPiPlot(igraph = igraph,method = "cluster_fast_greedy")
p <- res[[1]]
p
# ggsave("./cs2.pdf",p,width = 8,height = 7)扩展-关键OTU挑选
Hub节点是在网络中与其他节点连接较多的节点，Hub微生物就是与其他微生物联系较为紧密的微生物，可以称之为关键微生物（keystone）
hub = hub_score(igraph)$vector %>%
sort(decreasing = TRUE) %>%
head(5) %>%
as.data.frame()
colnames(hub) = "hub_sca"
ggplot(hub) +
geom_bar(aes(x = hub_sca,y = reorder(row.names(hub),hub_sca)),stat = "identity",fill = "#4DAF4A")
# ggsave("./cs2.pdf",width = 5,height = 4)对应随机网络构建和网络参数比对
result = random_Net_compate(igraph = igraph, type = "gnm", step = 100, netName = layout)
p1 = result[[1]]
sum_net = result[[4]]
p1
head(sum_net)## KO Means SD
## num.edges 1.078000e+03 1.078000e+03 0
## num.pos.edges 7.570000e+02 0.000000e+00 0
## num.neg.edges 3.210000e+02 0.000000e+00 0
## num.vertices 3.350000e+02 3.350000e+02 0
## connectance 1.926892e-02 1.926892e-02 0
## average.degree 6.435821e+00 6.435821e+00 0# ggsave("./cs3.pdf",p1,width = 5,height = 4)微生物网络流程
微生物组小网络：model_Gephi.2
使用network函数运行微生物网络全套分析：
- 使用OTU数量建议少于250个，如果OTU数量为250个，同时计算zipi，整个运算过程为3-5min。
data("ps16s")
path = "./result_micro_200/"
dir.create(path)
result = network(ps = ps16s,
N = 200,
layout_net = "model_Gephi.2",
r.threshold=0.6,
p.threshold=0.05,
label = FALSE,
path = path,
zipi = TRUE)## [1] "Group1"
## [1] "cor matrix culculating over"
## [1] "1"
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6## [1] "Group2"
## [1] "cor matrix culculating over"
## [1] "1"
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6## [1] "Group3"
## [1] "cor matrix culculating over"
## [1] "1"
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6# 多组网络绘制到一个面板
p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 48,height = 16,dpi = 72)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 48,height = 16)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)微生物大网络：model_maptree2
大网络运算时间会比较长，这里我没有计算zipi，用时5min完成全部运行。 N=0，代表用全部的OTU进行计算。
- 3000个OTU不计算zipi全套需要18min。
path = "./result_big_1000/"
dir.create(path)
result = network.2(ps = ps16s,
N = 1000,
big = TRUE,
maxnode = 5,
select_layout = TRUE,
layout_net = "model_maptree2",
r.threshold=0.4,
p.threshold=0.01,
label = FALSE,
path = path,
zipi = FALSE)## [1] "Group1"
## [1] "cor matrix culculating over"## [1] "Group2"
## [1] "cor matrix culculating over"## [1] "Group3"
## [1] "cor matrix culculating over"# 多组网络绘制到一个面板
p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
num= 3
# plotname1 = paste(path,"/network_all.jpg",sep = "")
# ggsave(plotname1, p,width = 16*num,height = 16,dpi = 72)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10*num,height = 10,limitsize = FALSE)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)微生物大网络：model_igraph
path = "./result_1000_igraph/"
dir.create(path)
# map = sample_data(ps)
# map$Group = "one"
# sample_data(ps16s) = map
result = network.2(ps = ps16s,
N = 1000,
big = TRUE,
maxnode = 5,
select_layout = TRUE,
layout_net = "model_igraph",
r.threshold=0.4,
p.threshold=0.01,
label = FALSE,
path = path,
zipi = FALSE)## [1] "Group1"
## [1] "cor matrix culculating over"## [1] "Group2"
## [1] "cor matrix culculating over"## [1] "Group3"
## [1] "cor matrix culculating over"# 多组网络绘制到一个面板
p = result[[1]]
# 全部样本网络参数比对
data = result[[2]]
num= 3
# plotname1 = paste(path,"/network_all.jpg",sep = "")
# ggsave(plotname1, p,width = 16*num,height = 16,dpi = 72)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 16*num,height = 16,limitsize = FALSE)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)微生物大网络：model_igraph2-network.i
model_igraph2在model_igraph的基础上去除了离散点，这对于重复数量少的研究比较有利，可以更加清楚的展示小样本网络。
network.i函数专门为model_igraph2算法而写，可以额外输出模块上色的网络。
path = "./result_1000_igraph2/"
dir.create(path)
# map = sample_data(ps)
# map$Group = "one"
# sample_data(ps16s) = map
result = network.i(ps = ps16s,
N = 1000,
r.threshold=0.9,
big = T,
select_layout = T,
method = "pearson",
scale = FALSE,
layout_net = "model_igraph2",
p.threshold=0.05,
label = FALSE,
path =path ,
zipi = F,
order = NULL
)## [1] "Group1"
## [1] "cor matrix culculating over"## [1] "Group2"
## [1] "cor matrix culculating over"## [1] "Group3"
## [1] "cor matrix culculating over"p1 = result[[1]]
p1
dat = result[[2]]
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(dat,tablename)
p = result[[5]]
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p1,width = 10,height = 3,limitsize = FALSE)
plotname1 = paste(path,"/network_all2.pdf",sep = "")
ggsave(plotname1, p,width = 25,height = 8,limitsize = FALSE)跨域网络流程
细菌，真菌，环境因子三者跨域网络
#--细菌和环境因子的网络#--------
Envnetplot<- paste("./16s_ITS_Env_network",sep = "")
dir.create(Envnetplot)
ps16s =ps16s%>% ggClusterNet::scale_micro()
psITS = psITS%>% ggClusterNet::scale_micro()
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
psITS= psITS,
NITS = 200,
N16s = 200)
map = phyloseq::sample_data(ps.merge)
head(map)## ID Group
## sample1 sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
#--环境因子导入
data1 = env
envRDA = data.frame(row.names = env$ID,env[,-1])
envRDA.s = vegan::decostand(envRDA,"hellinger")
data1[,-1] = envRDA.s
Gru = data.frame(ID = colnames(env)[-1],group = "env" )
head(Gru)## ID group
## 1 pH env
## 2 SOC env
## 3 TN env
## 4 NH4.N env
## 5 NO3.N env
## 6 AP envlibrary(sna)
library(ggClusterNet)
library(igraph)
result <- ggClusterNet::corBionetwork(ps = ps.merge,
N = 0,
r.threshold = 0.4, # 相关阈值
p.threshold = 0.05,
big = T,
group = "Group",
env = data1, # 环境指标表格
envGroup = Gru,# 环境因子分组文件表格
# layout = "fruchtermanreingold",
path = Envnetplot,# 结果文件存储路径
fill = "Phylum", # 出图点填充颜色用什么值
size = "igraph.degree", # 出图点大小用什么数据
scale = TRUE, # 是否要进行相对丰度标准化
bio = TRUE, # 是否做二分网络
zipi = F, # 是否计算ZIPI
step = 100, # 随机网络抽样的次数
width = 18,
label = TRUE,
height = 10
)## [1] "one"
## num [1:34, 1:18] 0.0474 0.0633 0.0188 0.0588 0.069 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:34] "pH" "SOC" "TN" "NH4.N" ...
## ..$ : chr [1:18] "sample1" "sample10" "sample11" "sample12" ...## [1] "1"## [1] "2"
## [1] "3"p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 15,height = 12,dpi = 72)
# plotname1 = paste(Envnetplot,"/network_all.png",sep = "")
# ggsave(plotname1, p,width = 10,height = 8,dpi = 72)
plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 15,height = 12)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)细菌-环境因子跨域网络
#--细菌-环境因子网络#-----
Envnetplot<- paste("./16s_Env_network",sep = "")
dir.create(Envnetplot)
# ps16s = ps0 %>% ggClusterNet::scale_micro()
psITS = NULL
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
psITS = NULL,
N16s = 500)
ps.merge## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 500 taxa and 18 samples ]
## sample_data() Sample Data: [ 18 samples by 2 sample variables ]
## tax_table() Taxonomy Table: [ 500 taxa by 8 taxonomic ranks ]map = phyloseq::sample_data(ps.merge)
head(map)## ID Group
## sample1 sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
data(env)
data1 = env
envRDA = data.frame(row.names = env$ID,env[,-1])
envRDA.s = vegan::decostand(envRDA,"hellinger")
data1[,-1] = envRDA.s
Gru = data.frame(ID = colnames(env)[-1],group = "env" )
head(Gru)## ID group
## 1 pH env
## 2 SOC env
## 3 TN env
## 4 NH4.N env
## 5 NO3.N env
## 6 AP env# library(sna)
# library(ggClusterNet)
# library(igraph)
result <- ggClusterNet::corBionetwork(ps = ps.merge,
N = 0,
r.threshold = 0.4, # 相关阈值
p.threshold = 0.05,
big = T,
group = "Group",
env = data1, # 环境指标表格
envGroup = Gru,# 环境因子分组文件表格
# layout = "fruchtermanreingold",
path = Envnetplot,# 结果文件存储路径
fill = "Phylum", # 出图点填充颜色用什么值
size = "igraph.degree", # 出图点大小用什么数据
scale = TRUE, # 是否要进行相对丰度标准化
bio = TRUE, # 是否做二分网络
zipi = F, # 是否计算ZIPI
step = 100, # 随机网络抽样的次数
width = 18,
label = TRUE,
height = 10
)## [1] "one"
## num [1:34, 1:18] 0.0474 0.0633 0.0188 0.0588 0.069 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:34] "pH" "SOC" "TN" "NH4.N" ...
## ..$ : chr [1:18] "sample1" "sample10" "sample11" "sample12" ...## [1] "1"## [1] "2"
## [1] "3"p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 13,height = 12,dpi = 72)
plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 13,height = 12)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)细菌和真菌跨域网络
#---细菌和真菌OTU网络-域网络-二分网络#-------
# 仅仅关注细菌和真菌之间的相关，不关注细菌内部和真菌内部相关
Envnetplot<- paste("./16S_ITS_network",sep = "")
dir.create(Envnetplot)
data(psITS)
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
psITS = psITS,
N16s = 500,
NITS = 500
)
ps.merge## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1000 taxa and 18 samples ]
## sample_data() Sample Data: [ 18 samples by 2 sample variables ]
## tax_table() Taxonomy Table: [ 1000 taxa by 8 taxonomic ranks ]map = phyloseq::sample_data(ps.merge)
head(map)## ID Group
## sample1 sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
data = data.frame(ID = c("fun_ASV_205","fun_ASV_316","fun_ASV_118"),
c("Coremode","Coremode","Coremode"))
# source("F:/Shared_Folder/Function_local/R_function/my_R_packages/ggClusterNet/R/corBionetwork2.R")
result <- corBionetwork(ps = ps.merge,
N = 0,
lab = data,
r.threshold = 0.8, # 相关阈值
p.threshold = 0.05,
group = "Group",
# env = data1, # 环境指标表格
# envGroup = Gru,# 环境因子分组文件表格
# layout = "fruchtermanreingold",
path = Envnetplot,# 结果文件存储路径
fill = "Phylum", # 出图点填充颜色用什么值
size = "igraph.degree", # 出图点大小用什么数据
scale = TRUE, # 是否要进行相对丰度标准化
bio = TRUE, # 是否做二分网络
zipi = F, # 是否计算ZIPI
step = 100, # 随机网络抽样的次数
width = 12,
label = TRUE,
height = 10,
big = TRUE,
select_layout = TRUE,
layout_net = "model_maptree2",
clu_method = "cluster_fast_greedy"
)## [1] "one"## [1] "1"## [1] "2"
## [1] "3"tem <- model_maptree(cor =result[[5]],
method = "cluster_fast_greedy",
seed = 12
)
node_model = tem[[2]]
head(node_model)## ID group degree
## 1 fun_ASV_246 40 35
## 2 fun_ASV_770 40 33
## 3 fun_ASV_939 40 30
## 4 fun_ASV_233 41 22
## 5 fun_ASV_76 43 22
## 6 fun_ASV_268 19 19p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 20,height = 19)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)
tablename <- paste(Envnetplot,"/node_model_imformation",".csv",sep = "")
write.csv(node_model,tablename)
tablename <- paste(Envnetplot,"/nodeG_plot",".csv",sep = "")
write.csv(result[[4]],tablename)
tablename <- paste(Envnetplot,"/edge_plot",".csv",sep = "")
write.csv(result[[3]],tablename)
tablename <- paste(Envnetplot,"/cor_matrix",".csv",sep = "")
write.csv(result[[5]],tablename)细菌真菌任意分类水平跨域网络
#--细菌真菌任意分类水平跨域网络#----------
library(tidyverse)
# res1path = "result_and_plot/Micro_and_other_index_16s/"
Envnetplot<- paste("./16S_ITS_network_Genus",sep = "")
dir.create(Envnetplot)
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- merge16S_ITS(ps16s = ps16s,
psITS = psITS,
N16s = 500,
NITS = 500
)
ps.merge## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1000 taxa and 18 samples ]
## sample_data() Sample Data: [ 18 samples by 2 sample variables ]
## tax_table() Taxonomy Table: [ 1000 taxa by 8 taxonomic ranks ]map = phyloseq::sample_data(ps.merge)
head(map)## ID Group
## sample1 sample1 Group1
## sample10 sample10 Group2
## sample11 sample11 Group2
## sample12 sample12 Group2
## sample13 sample13 Group3
## sample14 sample14 Group3map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
tem.0 = ps.merge %>% tax_glom_wt(ranks = "Genus")
tax = tem.0 %>% vegan_tax() %>%
as.data.frame()
head(tax)## filed Kingdom Phylum Class
## Acaulium fun Fungi Ascomycota Sordariomycetes
## Acidocella bac Bacteria Proteobacteria Alphaproteobacteria
## Acrobeloides fun Metazoa Nematoda Chromadorea
## Actinomadura bac Bacteria Actinobacteria Actinobacteria
## Altererythrobacter bac Bacteria Proteobacteria Alphaproteobacteria
## Alternaria fun Fungi Ascomycota Dothideomycetes
## Order Family Genus
## Acaulium Microascales Microascaceae Acaulium
## Acidocella Rhodospirillales Acetobacteraceae Acidocella
## Acrobeloides Rhabditida Cephalobidae Acrobeloides
## Actinomadura Actinomycetales Thermomonosporaceae Actinomadura
## Altererythrobacter Sphingomonadales Erythrobacteraceae Altererythrobacter
## Alternaria Pleosporales Pleosporaceae Alternariadata = data.frame(ID = c("Acaulium","Acidocella","Acrobeloides"),
c("Acaulium","Acidocella","Acrobeloides"))
result <- corBionetwork(ps = tem.0,
N = 0,
lab = data,
r.threshold = 0.6, # 相关阈值
p.threshold = 0.05,
group = "Group",
# env = data1, # 环境指标表格
# envGroup = Gru,# 环境因子分组文件表格
# layout = "fruchtermanreingold",
path = Envnetplot,# 结果文件存储路径
fill = "Phylum", # 出图点填充颜色用什么值
size = "igraph.degree", # 出图点大小用什么数据
scale = TRUE, # 是否要进行相对丰度标准化
bio = TRUE, # 是否做二分网络
zipi = F, # 是否计算ZIPI
step = 100, # 随机网络抽样的次数
width = 12,
label = TRUE,
height = 10,
big = TRUE,
select_layout = TRUE,
layout_net = "model_maptree2",
clu_method = "cluster_fast_greedy"
)## [1] "one"## [1] "1"## [1] "2"
## [1] "3"tem <- model_maptree(cor =result[[5]],
method = "cluster_fast_greedy",
seed = 12
)
node_model = tem[[2]]
head(node_model)## ID group degree
## 1 Microascus 12 13
## 2 Scopulariopsis 13 12
## 3 Arcopilus 10 10
## 4 Cladosporium 13 9
## 5 Conocybe 13 8
## 6 Hydrogonium 6 8library(WGCNA)
otu = tem.0 %>% vegan_otu() %>%
as.data.frame()
node_model = node_model[match(colnames(otu),node_model$ID),]
MEList = moduleEigengenes(otu, colors = node_model$group)
MEs = MEList$eigengenes
nGenes = ncol(otu)
nSamples = nrow(otu)
moduleTraitCor = cor(MEs, envRDA, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
#sizeGrWindow(10,6)
pdf(file=paste(Envnetplot,"/","Module-env_relationships.pdf",sep = ""),width=10,height=6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(envRDA),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
dev.off()## png
## 2p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10,height = 10)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)
tablename <- paste(Envnetplot,"/nodeG_plot",".csv",sep = "")
write.csv(result[[4]],tablename)
tablename <- paste(Envnetplot,"/edge_plot",".csv",sep = "")
write.csv(result[[3]],tablename)
tablename <- paste(Envnetplot,"/cor_matrix",".csv",sep = "")
write.csv(result[[5]],tablename)网络稳定性（抗干扰性）
模块相似度
#--网络稳定性:模块相似性------
data(ps)
library(tidyfst)
res = module.compare.m(
ps = ps,
Top = 500,
degree = TRUE,
zipi = FALSE,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman",
padj = F,
n = 3)## [1] 80 8
## [1] 159 8
## [1] 224 8
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 0
## [1] 1
## [1] 2
## [1] 3
## [1] 4#不同分组使用一个圆圈展示，圆圈内一个点代表一个模块，相连接的模块代表了相似的模块。
p = res[[1]]
p
#--提取模块的OTU，分组等的对应信息
dat = res[[2]]
head(dat)## ID group Group
## 1 ASV_256 KOmodel_1 KO
## 2 ASV_142 KOmodel_1 KO
## 3 ASV_362 KOmodel_1 KO
## 4 ASV_89 KOmodel_1 KO
## 5 ASV_419 KOmodel_1 KO
## 6 ASV_352 KOmodel_1 KO#模块相似度结果表格
dat2 = res[[3]]
head(dat2)## module1 module2 Both P1A2 P2A1 A1A2
## KOmodel_41-OEmodel_2 KOmodel_41 OEmodel_2 1 2 6 500
## KOmodel_62-OEmodel_1 KOmodel_62 OEmodel_1 1 2 6 500
## KOmodel_34-OEmodel_1 KOmodel_34 OEmodel_1 1 2 6 500
## KOmodel_22-OEmodel_4 KOmodel_22 OEmodel_4 1 3 5 500
## KOmodel_26-OEmodel_8 KOmodel_26 OEmodel_8 1 3 4 501
## KOmodel_17-OEmodel_10 KOmodel_17 OEmodel_10 1 4 4 500
## p_raw p_adj
## KOmodel_41-OEmodel_2 0.0407716775257314 0.0407716775257314
## KOmodel_62-OEmodel_1 0.0407716775257314 0.0407716775257314
## KOmodel_34-OEmodel_1 0.0407716775257314 0.0407716775257314
## KOmodel_22-OEmodel_4 0.0464588019672787 0.0464588019672787
## KOmodel_26-OEmodel_8 0.0388304723830171 0.0388304723830171
## KOmodel_17-OEmodel_10 0.0483470023594227 0.0483470023594227网络鲁棒性（随机去除节点）
这里通过随机去除部分OTU，计算网络鲁棒性，代表网络抗干扰的能力。 按照0.05的步长，每次去除5%的文生物，重新计算鲁棒性，知道最终全部去除。 如果有分组列Group，则会按照分组进行鲁棒性计算，并且区分颜色绘制到一个面板中。 计算鲁棒性这里使用丰度加成权重和不加权两种方式，左边是不加权，后侧是加权的结果。 这里步长不可以修改，因为修改好像也没什么意思。
#--随即取出任意比例节点-网络鲁棒性#---------
res = Robustness.Random.removal(ps = ps,
Top = 500,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman"
)
p = res[[1]]
p
#提取数据
dat = res[[2]]
head(dat)## Proportion.removed remain.mean remain.sd remain.se weighted Group
## 1 0.05 0.5408197 0.01310521 0.001310521 weighted KO
## 2 0.10 0.5030601 0.01536428 0.001536428 weighted KO
## 3 0.15 0.4674044 0.01755406 0.001755406 weighted KO
## 4 0.20 0.4325683 0.01966370 0.001966370 weighted KO
## 5 0.25 0.3944809 0.02034929 0.002034929 weighted KO
## 6 0.30 0.3542623 0.02243074 0.002243074 weighted KO网络鲁棒性（去除关键节点）
#---去除关键节点-网络鲁棒性#------
res= Robustness.Targeted.removal(ps = ps,
Top = 500,
degree = TRUE,
zipi = FALSE,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman")
p = res[[1]]
p
#提取数据
dat = res[[2]]负相关比例
#--计算负相关的比例#----
res = negative.correlation.ratio(ps = ps,
Top = 500,
degree = TRUE,
zipi = FALSE,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman")
p = res[[1]]
p
dat = res[[2]]
#-负相关比例数据
head(dat)## ID ratio
## 1 KO 41.15139
## 2 OE 34.19023
## 3 WT 39.24051# dir.create("./negative_correlation_ratio/")
# write.csv(dat,
# paste("./negative_correlation_ratio/","_negative_ratio_network.csv",sep = ""))组成稳定性
treat = ps %>% sample_data()
treat$pair = paste( "A",c(rep(1:6,3)),sep = "")
head(treat)## Group Date Site Sequencing Platform Species Batch
## KO1 KO 2017/6/30 Beijing BGI HiSeq2500 Arabidopsis 1
## KO2 KO 2017/6/30 Beijing BGI HiSeq2500 Arabidopsis 1
## KO3 KO 2017/7/2 Sanya BGI HiSeq2500 Arabidopsis 1
## KO4 KO 2017/7/2 Sanya BGI HiSeq2500 Arabidopsis 1
## KO5 KO 2017/7/4 Harbin BGI HiSeq2500 Arabidopsis 1
## KO6 KO 2017/7/5 Harbin BGI HiSeq2500 Arabidopsis 1
## BarcodeSequence LinkerPrimerSequence ReversePrimer ID pair
## KO1 ACGCTCGACA AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO1 A1
## KO2 ATCAGACACG AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO2 A2
## KO3 ATATCGCGAG AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO3 A3
## KO4 CACGAGACAG AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO4 A4
## KO5 CTCGCGTGTC AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO5 A5
## KO6 TAGTATCAGC AACMGGATTAGATACCCKG ACGTCATCCCCACCTTCC KO6 A6sample_data(ps) = treat
#一般性的处理，没有时间梯度的，这里设定time为F，意味着每两个群落的组合都进行比较
res = community.stability( ps = ps,time = FALSE)
p = res[[1]]
p
dat = res[[2]]
# 如果是时间序列，会按照时间的单一流动性方向进行比较，自然就不是两两比对了。
res = community.stability( ps = ps,time = TRUE)
p = res[[1]]
p
网络抗毁性
网络抗毁性:使用了自然连通度的概念，用于反映网络稳定性.具体而言，通过在构建好的网络中移除里面的节点，并计算自然连通度。这样随着取出点的数量逐渐增多，就可以计算得到一连串的网络连通度。这个算法最先来自PNAS的文章：“Reduced microbial stability in the active layer is associated with carbon loss under alpine permafrost degradation”
library(tidyfst)
library("pulsar")
library(ggClusterNet)
library(phyloseq)
library(tidyverse)
res = natural.con.microp (
ps = ps,
Top = 200,
r.threshold= 0.5,
p.threshold=0.05,
method = "spearman",
norm = F,
end = 150,# 小于网络包含的节点数量
start = 0,
con.method = "pulsar"
)
p = res[[1]]
p
dat = res[[2]]网络模块化
展示网络模块
这里我们通过指定模块
library(ggClusterNet)
library(igraph)
library(ggClusterNet)
library(tidyverse)
library(phyloseq)
data(ps)
resu = module_display.2(
pst = ps,
r.threshold= 0.6,
p.threshold=0.05,
select.mod = c("model_1","model_2","model_3","model_4"),#选择指定模块可视化
Top = 500,
num = 5, # 模块包含OTU数量少于5个的不展示,
leg.col = 3
)
# 全部模块输出展示
p1 = resu[[1]]
p1
#--提取率相关矩阵
dat0 = resu[[4]]
#--提取模块信息：这里区分好不同模块的OTU，以及每个模块包含的网络边的数量
dat = resu[[5]]
head(dat)## ID group degree
## 1 ASV_66 model_3 46
## 2 ASV_294 model_3 40
## 3 ASV_326 model_3 34
## 4 ASV_60 model_2 31
## 5 ASV_2 model_3 31
## 6 ASV_163 model_3 30#--提取网络的边和节点表格
dat2 = resu[[6]]
names(dat2)## [1] "edge" "node"指定的目标模块展示
p2 = resu[[2]]
p2
#按照模块着色，模块之间的边着色为灰色展示整个网络。
p2 = resu[[3]]
p2
#--注意这个报错信息，为展示的空间不够大，拉一拉展示框即可
# Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto), :
# Viewport has zero dimension(s)基于模块的微生物组成分析
对微生物数据分好模块，使用模块分类的数据导入module_composition函数，即可作指定模块的物种组成。
#--这里我们指定三个模块，绘制这三个模块物种组成图表
select.mod = c("model_1","model_2","model_3")
#--模块信息，共有三列，第一列是OTU，第二列是吗模块的分类，第三列是模块捏的边数量
mod1 = resu$mod.groups %>% filter(group %in% select.mod)
head(mod1)## ID group degree
## 1 ASV_66 model_3 46
## 2 ASV_294 model_3 40
## 3 ASV_326 model_3 34
## 4 ASV_60 model_2 31
## 5 ASV_2 model_3 31
## 6 ASV_163 model_3 30#-计算模块物种组成信息
pst = ps %>%
filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
scale_micro("rela") %>%
# subset_samples.wt("Group" ,id3[m]) %>%
filter_OTU_ps(500)
#-选择模块的微生物组成，通过参数j选择不同的分类水平。
res = module_composition(pst = pst,mod1 = mod1,j = "Family")
p3 = res[[1]]
p3
#按照属水平进行绘制
res = module_composition(pst = pst,mod1 = mod1,j = "Genus")
p3 = res[[1]]
p3
#---提取出图数据物种组成数据导出
ps.t = res[[3]]
otu = ps.t %>% vegan_otu() %>% t() %>%
as.data.frame()
# head(otu)
tax = ps.t %>% vegan_tax() %>%
as.data.frame()
# head(tax)
tab.res = cbind(otu,tax)
# head(tab.res)
# 没有标准化的出图数据
tab.res.2 = res[[4]]$bundance
# head(tab.res.2)
# 标准化到100%的出图数据
tab.res.3 = res[[4]]$relaabundance
# head(tab.res.3)模块计算alpha多样性
这里将指定模块和ps对象输入函数module_alpha即可计算三种alpha多样性，这里调用了EasyStat包，使用非参数检验检测了alpha多样性的显著性。
注意：这里的ps对象需要原始count数据，不可用标准化后的数据。
#--这里我们指定三个模块，绘制这三个模块物种组成图表
select.mod = c("model_1","model_2","model_3")
#--模块信息，两列，第一列是OTU，第二列是模块的分类
mod1 = resu$mod.groups %>%
dplyr::select(ID,group) %>%
dplyr::filter(group %in% select.mod)
head(mod1)## ID group
## 1 ASV_66 model_3
## 2 ASV_294 model_3
## 3 ASV_326 model_3
## 4 ASV_60 model_2
## 5 ASV_2 model_3
## 6 ASV_163 model_3#-选择模块的多样性
result = module_alpha (
ps = ps,
mod1 = mod1)
p5 = result[[1]]
p5
#--alpha多样性数据
plotd = result[[4]]$alpha
# alpha多样性显著性表格
sigtab = result[[4]]$sigtab基于单个样本计算网络属性
这里我们需要输入微生物组数据的ps对象，其次，需要输入这组微生物的相关矩阵；然后netproperties.sample函数会按照单个样本计算网络属性；这里一共有15中网络属性。
这里计算出来后便可以和其他指标进行相关性分析；
library(ggClusterNet)
library(igraph)
library(ggClusterNet)
library(tidyverse)
library(phyloseq)
data(ps)
#-计算模块物种组成信息
pst = ps %>%
filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
scale_micro("rela") %>%
# subset_samples.wt("Group" ,id3[m]) %>%
filter_OTU_ps(500)
# 选择和网络属性做相关的指标
select.env = "env1"
#--内置数据无需导入
env = env1 %>%
as.data.frame() %>%
rownames_to_column("id")
dat.f = netproperties.sample(pst = pst,cor = cor)
# head(dat.f)
dat.f$id = row.names(dat.f)
dat.f = dat.f %>% dplyr:: select(id,everything())
subenv = env %>% dplyr::select(id,everything()) %>% dplyr::select(id,select.env )
tab = dat.f %>% left_join(subenv,by = "id")
mtcars2 = reshape2::melt(tab, id.vars=c(select.env,"id"))
lab = mean(mtcars2[,select.env])
preoptab = tab
p0_1 = ggplot2::ggplot(mtcars2,aes(x= value,!!sym(select.env), colour=variable)) +
ggplot2::geom_point() +
ggpubr::stat_cor(label.y=lab*1.1)+
ggpubr::stat_regline_equation(label.y=lab*1.1,vjust = 2) +
facet_wrap(~variable, scales="free_x") +
geom_smooth(aes(value,!!sym(select.env), colour=variable), method=lm, se=T)+
theme_classic()
p0_1
计算模块的Zscore值
按照样本计算模块的ZS值，方便施用这一指标和其他相关指标进行相关性分析。
dat = ZS.net.module(
pst = ps,
Top = 500,
r.threshold= 0.8,
p.threshold=0.05,
select.mod = NULL
)
# head(dat)
map =sample_data(ps)
map$id = row.names(map)
map = map[,c("id","Group")]
data = map %>%
as.tibble() %>%
inner_join(dat,by = "id") %>%
dplyr::rename(group = Group)
colnames(env)[1] = "id"
subenv = env %>% dplyr::select(id,everything()) %>% dplyr::select(id,select.env )
# head(data)
tab = data %>% left_join(subenv,by = "id")
modenv = tab
mtcars2 = reshape2::melt(tab, id.vars=c(select.env,"group","id"))
mtcars2$variable## [1] model_1 model_1 model_1 model_1 model_1 model_1 model_1
## [8] model_1 model_1 model_1 model_1 model_1 model_1 model_1
## [15] model_1 model_1 model_1 model_1 model_3 model_3 model_3
## [22] model_3 model_3 model_3 model_3 model_3 model_3 model_3
## [29] model_3 model_3 model_3 model_3 model_3 model_3 model_3
## [36] model_3 model_2 model_2 model_2 model_2 model_2 model_2
## [43] model_2 model_2 model_2 model_2 model_2 model_2 model_2
## [50] model_2 model_2 model_2 model_2 model_2 model_4 model_4
## [57] model_4 model_4 model_4 model_4 model_4 model_4 model_4
## [64] model_4 model_4 model_4 model_4 model_4 model_4 model_4
## [71] model_4 model_4 model_7 model_7 model_7 model_7 model_7
## [78] model_7 model_7 model_7 model_7 model_7 model_7 model_7
## [85] model_7 model_7 model_7 model_7 model_7 model_7 model_31
## [92] model_31 model_31 model_31 model_31 model_31 model_31 model_31
## [99] model_31 model_31 model_31 model_31 model_31 model_31 model_31
## [106] model_31 model_31 model_31 model_12 model_12 model_12 model_12
## [113] model_12 model_12 model_12 model_12 model_12 model_12 model_12
## [120] model_12 model_12 model_12 model_12 model_12 model_12 model_12
## [127] model_5 model_5 model_5 model_5 model_5 model_5 model_5
## [134] model_5 model_5 model_5 model_5 model_5 model_5 model_5
## [141] model_5 model_5 model_5 model_5 model_23 model_23 model_23
## [148] model_23 model_23 model_23 model_23 model_23 model_23 model_23
## [155] model_23 model_23 model_23 model_23 model_23 model_23 model_23
## [162] model_23 model_19 model_19 model_19 model_19 model_19 model_19
## [169] model_19 model_19 model_19 model_19 model_19 model_19 model_19
## [176] model_19 model_19 model_19 model_19 model_19 model_14 model_14
## [183] model_14 model_14 model_14 model_14 model_14 model_14 model_14
## [190] model_14 model_14 model_14 model_14 model_14 model_14 model_14
## [197] model_14 model_14 model_6 model_6 model_6 model_6 model_6
## [204] model_6 model_6 model_6 model_6 model_6 model_6 model_6
## [211] model_6 model_6 model_6 model_6 model_6 model_6 model_15
## [218] model_15 model_15 model_15 model_15 model_15 model_15 model_15
## [225] model_15 model_15 model_15 model_15 model_15 model_15 model_15
## [232] model_15 model_15 model_15 model_29 model_29 model_29 model_29
## [239] model_29 model_29 model_29 model_29 model_29 model_29 model_29
## [246] model_29 model_29 model_29 model_29 model_29 model_29 model_29
## [253] model_13 model_13 model_13 model_13 model_13 model_13 model_13
## [260] model_13 model_13 model_13 model_13 model_13 model_13 model_13
## [267] model_13 model_13 model_13 model_13 model_8 model_8 model_8
## [274] model_8 model_8 model_8 model_8 model_8 model_8 model_8
## [281] model_8 model_8 model_8 model_8 model_8 model_8 model_8
## [288] model_8 model_22 model_22 model_22 model_22 model_22 model_22
## [295] model_22 model_22 model_22 model_22 model_22 model_22 model_22
## [302] model_22 model_22 model_22 model_22 model_22 model_25 model_25
## [309] model_25 model_25 model_25 model_25 model_25 model_25 model_25
## [316] model_25 model_25 model_25 model_25 model_25 model_25 model_25
## [323] model_25 model_25 model_30 model_30 model_30 model_30 model_30
## [330] model_30 model_30 model_30 model_30 model_30 model_30 model_30
## [337] model_30 model_30 model_30 model_30 model_30 model_30 model_26
## [344] model_26 model_26 model_26 model_26 model_26 model_26 model_26
## [351] model_26 model_26 model_26 model_26 model_26 model_26 model_26
## [358] model_26 model_26 model_26 model_20 model_20 model_20 model_20
## [365] model_20 model_20 model_20 model_20 model_20 model_20 model_20
## [372] model_20 model_20 model_20 model_20 model_20 model_20 model_20
## [379] model_9 model_9 model_9 model_9 model_9 model_9 model_9
## [386] model_9 model_9 model_9 model_9 model_9 model_9 model_9
## [393] model_9 model_9 model_9 model_9 model_16 model_16 model_16
## [400] model_16 model_16 model_16 model_16 model_16 model_16 model_16
## [407] model_16 model_16 model_16 model_16 model_16 model_16 model_16
## [414] model_16 model_10 model_10 model_10 model_10 model_10 model_10
## [421] model_10 model_10 model_10 model_10 model_10 model_10 model_10
## [428] model_10 model_10 model_10 model_10 model_10 model_28 model_28
## [435] model_28 model_28 model_28 model_28 model_28 model_28 model_28
## [442] model_28 model_28 model_28 model_28 model_28 model_28 model_28
## [449] model_28 model_28 model_33 model_33 model_33 model_33 model_33
## [456] model_33 model_33 model_33 model_33 model_33 model_33 model_33
## [463] model_33 model_33 model_33 model_33 model_33 model_33 model_18
## [470] model_18 model_18 model_18 model_18 model_18 model_18 model_18
## [477] model_18 model_18 model_18 model_18 model_18 model_18 model_18
## [484] model_18 model_18 model_18 model_21 model_21 model_21 model_21
## [491] model_21 model_21 model_21 model_21 model_21 model_21 model_21
## [498] model_21 model_21 model_21 model_21 model_21 model_21 model_21
## [505] model_17 model_17 model_17 model_17 model_17 model_17 model_17
## [512] model_17 model_17 model_17 model_17 model_17 model_17 model_17
## [519] model_17 model_17 model_17 model_17 model_11 model_11 model_11
## [526] model_11 model_11 model_11 model_11 model_11 model_11 model_11
## [533] model_11 model_11 model_11 model_11 model_11 model_11 model_11
## [540] model_11 model_27 model_27 model_27 model_27 model_27 model_27
## [547] model_27 model_27 model_27 model_27 model_27 model_27 model_27
## [554] model_27 model_27 model_27 model_27 model_27 model_32 model_32
## [561] model_32 model_32 model_32 model_32 model_32 model_32 model_32
## [568] model_32 model_32 model_32 model_32 model_32 model_32 model_32
## [575] model_32 model_32 model_24 model_24 model_24 model_24 model_24
## [582] model_24 model_24 model_24 model_24 model_24 model_24 model_24
## [589] model_24 model_24 model_24 model_24 model_24 model_24 mother_no
## [596] mother_no mother_no mother_no mother_no mother_no mother_no mother_no
## [603] mother_no mother_no mother_no mother_no mother_no mother_no mother_no
## [610] mother_no mother_no mother_no
## 34 Levels: model_1 model_3 model_2 model_4 model_7 model_31 ... mother_nohead(mtcars2)## env1 group id variable value
## 1 7.67 KO KO1 model_1 7.1549454
## 2 7.61 KO KO2 model_1 -0.3726327
## 3 7.15 KO KO3 model_1 10.6994573
## 4 6.89 KO KO4 model_1 10.3945282
## 5 7.39 KO KO5 model_1 35.8009267
## 6 7.01 KO KO6 model_1 9.5297211lab = mean(mtcars2[,select.env])
p1_1 = ggplot2::ggplot(mtcars2,aes(x= value,!!sym(select.env), colour=variable)) +
ggplot2::geom_point() +
ggpubr::stat_cor(label.y=lab*1.1)+
ggpubr::stat_regline_equation(label.y=lab*1.1,vjust = 2) +
facet_wrap(~variable, scales="free_x") +
geom_smooth(aes(value,!!sym(select.env), colour=variable), method=lm, se=T)+
theme_classic()
p1_1
模块ZS值，网络属性，环境因子相关联合到一起
#网络模块和环境因子相关-网络指标和环境相关
dat.1 = env1 %>% rownames_to_column("id")
res = net.property.module.env (
pst = pst,
Top = 500,
r.threshold= 0.8,
p.threshold=0.05,
env = dat.1,
select.mod = c("model_1","model_2","model_3"),
select.env = "env1")
#--模块和环境因子之间相关
p9 = res[[1]]
p9
#-网络属性和环境因子相关
p9 = res[[2]]
p9
时空组:网络相关分析
时空组：多组网络展示-网络性质稳定性等
这部分更新与2023年1月27日完成，可能会存在bug等，如有问题即使告知我解决；
这部分更新主要是为了解决一下几个问题： - 多个网络在一张图上需要调整映射颜色相同 - 分组并不是只有一列，很多时候有时间序列，有空间部位等，优化参数使得时间序列等多个分组的研究展示更加顺畅；
rm(list=ls())
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)
library(tidyfst)
ps.st = readRDS("./ps_TS.rds")
ps.st## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2432 taxa and 108 samples ]
## sample_data() Sample Data: [ 108 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 2432 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
## refseq() DNAStringSet: [ 2432 reference sequences ]#时空组网络-分面网络图-解决填充颜色不一致问题#----
res = Facet.network (
ps.st= ps.st,# phyloseq对象
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time",# 分组3
ord.g1 = c("WT","KO","OE"),# 排序顺序
ord.g2 = c("B","R") ,# 排序顺序
ord.g3 = c("T1","T2","T3") ,# 排序顺序
order = "time", # 出图每行代表的变量
fill = "Phylum",
size = "igraph.degree",
layout_net = "model_maptree2",
r.threshold=0.8,
p.threshold=0.01,
method = "spearman",
select_layout = TRUE,
clu_method = "cluster_fast_greedy",
maxnode = 5
)## [1] "B" "R"
## [1] "T1" "T2" "T3"
## [1] "WT" "KO" "OE"p = res[[1]]
p
# ggsave("cs1.pdf",p,width = 5*5,height = 4*3)
#--时空组网络-网络性质#-------
dat = net.pro.ts(dat = res[[2]][[1]])
head(dat)## B.T1.WT B.T1.KO
## num.edges(L) "81" "96"
## num.pos.edges "41" "60"
## num.neg.edges "40" "36"
## num.vertices(n) "100" "106"
## Connectance(edge_density) "0.0163636363636364" "0.0172506738544474"
## average.degree(Average K) "1.62" "1.81132075471698"
## B.T1.OE B.T2.WT
## num.edges(L) "66" "65"
## num.pos.edges "46" "35"
## num.neg.edges "20" "30"
## num.vertices(n) "96" "91"
## Connectance(edge_density) "0.0144736842105263" "0.0158730158730159"
## average.degree(Average K) "1.375" "1.42857142857143"
## B.T2.KO B.T2.OE
## num.edges(L) "106" "78"
## num.pos.edges "67" "55"
## num.neg.edges "39" "23"
## num.vertices(n) "105" "109"
## Connectance(edge_density) "0.0194139194139194" "0.0132517838939857"
## average.degree(Average K) "2.01904761904762" "1.43119266055046"
## B.T3.WT B.T3.KO
## num.edges(L) "73" "92"
## num.pos.edges "37" "55"
## num.neg.edges "36" "37"
## num.vertices(n) "94" "103"
## Connectance(edge_density) "0.0167009837565774" "0.0175138016371597"
## average.degree(Average K) "1.5531914893617" "1.78640776699029"
## B.T3.OE R.T1.WT
## num.edges(L) "77" "72"
## num.pos.edges "53" "36"
## num.neg.edges "24" "36"
## num.vertices(n) "105" "96"
## Connectance(edge_density) "0.0141025641025641" "0.0157894736842105"
## average.degree(Average K) "1.46666666666667" "1.5"
## R.T1.KO R.T1.OE
## num.edges(L) "95" "82"
## num.pos.edges "51" "63"
## num.neg.edges "44" "19"
## num.vertices(n) "107" "108"
## Connectance(edge_density) "0.01675189560924" "0.0141917618553133"
## average.degree(Average K) "1.77570093457944" "1.51851851851852"
## R.T2.WT R.T2.KO
## num.edges(L) "67" "107"
## num.pos.edges "33" "60"
## num.neg.edges "34" "47"
## num.vertices(n) "95" "106"
## Connectance(edge_density) "0.0150055991041433" "0.0192273135669362"
## average.degree(Average K) "1.41052631578947" "2.0188679245283"
## R.T2.OE R.T3.WT
## num.edges(L) "81" "72"
## num.pos.edges "60" "39"
## num.neg.edges "21" "33"
## num.vertices(n) "104" "97"
## Connectance(edge_density) "0.0151232262882748" "0.0154639175257732"
## average.degree(Average K) "1.55769230769231" "1.48453608247423"
## R.T3.KO R.T3.OE
## num.edges(L) "94" "80"
## num.pos.edges "56" "57"
## num.neg.edges "38" "23"
## num.vertices(n) "99" "108"
## Connectance(edge_density) "0.0193774479488765" "0.0138456213222568"
## average.degree(Average K) "1.8989898989899" "1.48148148148148"#--时空组的网络稳定性
#--时空组-稳定性1模块比对#--------
res = module.compare.m.ts(ps.st = ps.st, N = 200,
degree = TRUE,zipi = FALSE,r.threshold= 0.8,
p.threshold=0.05,method = "spearman",
padj = F,n = 3,
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time",# 分组3
zoom = 0.3,# 控制小圈大小
b.x = 1.3,
b.y = 1.3)## [1] 1 8
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 2 8
## [1] 2 8
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 2 8
## NULL
## NULL
## NULL
## [1] 2 8
## [1] 2 8
## NULL
## [1] 1 8
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## [1] 2 8
## NULL
## [1] 1 8
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 1 8
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 2 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 1 8
## [1] 1 8
## [1] 2 8
## NULL
## [1] 2 8
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## NULL
## [1] 1 8
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 0
## [1] 1
## [1] 2
## [1] 0
## [1] 1
## [1] 2p = res[[1]]
p
#提取数据
dat = res[[2]]
#-作图数据提取
dat = res[[3]]
#时空组-稳定性2 鲁棒性随即取出#-----
res = Robustness.Random.removal.ts(
ps.st = ps.st,
N = 200,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman",
order = "time",
g1 = "Group",
g2 = "space",
g3 = "time")
#
res[[1]][[1]]
res[[1]][[2]]
#时空组-稳定性2-2 鲁棒性随即取出#-----
res = Robustness.Targeted.removal.ts(
ps.st,
N = 200,
degree = TRUE,
zipi = FALSE,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman",
order = "time",
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time")
res[[1]][[1]]
res[[1]][[2]]
#时空组-稳定性3负相关比例#-------
res = negative.correlation.ratio.ts(
ps.st,
N = 200,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman",
order = "time",
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time",# 分组3
ord.g1 =NULL,# 排序顺序
ord.g2 = NULL, # 排序顺序
ord.g3 = NULL# 排序顺序
)
# 导出图片
res[[1]]
# 导出数据
dat = res[[2]]
head(dat)## ID ratio group time space label Group
## 1 R.T1.KO 46.31579 KO T1 R T1.R R.T1.KO
## 2 R.T1.OE 23.17073 OE T1 R T1.R R.T1.OE
## 3 R.T1.WT 50.00000 WT T1 R T1.R R.T1.WT
## 4 R.T2.KO 43.92523 KO T2 R T2.R R.T2.KO
## 5 R.T2.OE 25.92593 OE T2 R T2.R R.T2.OE
## 6 R.T2.WT 50.74627 WT T2 R T2.R R.T2.WT#时空组-稳定性4 网络易损性#------
library("pulsar")
res = natural.con.microp.ts(
ps.st,
N = 200,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman",
order = "time",
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time",# 分组3
ord.g1 =NULL,# 排序顺序
ord.g2 = NULL, # 排序顺序
ord.g3 = NULL,# 排序顺序
norm = F,
end = N -2,
start = 0,
con.method = "pulsar"
)
res[[1]]
#时空组-稳定性5 组成稳定性#-------
res = community.stability.ts(
ps.st = ps.st,
N = 200,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman",
order = "space",
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time",# 分组3
map.art = NULL, # 人工输入的分组 默认为NULL
time = F,# 稳定性是否有时间序列
ord.map = TRUE# map文件是否是已经按照pair要求进行了排序
)
res[[1]]
# res[[2]]时空组:多网络单独绘制却填充相同颜色
#---时空组双网络手动填充相同颜色#-------
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)
ps.st = readRDS("./ps_TS.rds")
ps.st## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2432 taxa and 108 samples ]
## sample_data() Sample Data: [ 108 samples by 4 sample variables ]
## tax_table() Taxonomy Table: [ 2432 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
## refseq() DNAStringSet: [ 2432 reference sequences ]#-----第一组网络绘制
res = Facet.network(
ps.st= ps.st,# phyloseq对象
g1 = "Group",# 分组1
g2 = "space",# 分组2
g3 = "time",# 分组3
ord.g1 = c("WT","KO","OE"),# 排序顺序
ord.g2 = c("B","R") ,# 排序顺序
ord.g3 = c("T1","T2","T3") ,# 排序顺序
order = "time", # 出图每行代表的变量
fill = "Phylum",
size = "igraph.degree",
layout_net = "model_maptree2",
r.threshold=0.8,
p.threshold=0.01,
method = "spearman",
select_layout = TRUE,
clu_method = "cluster_fast_greedy",
maxnode = 5
)## [1] "B" "R"
## [1] "T1" "T2" "T3"
## [1] "WT" "KO" "OE"p = res[[1]]
p
#--指定颜色映射，为了多个图使用同一个颜色#-------
#--设定的参数一致
fill = "Phylum"
size = "igraph.degree"
maxnode = 5
row.num = 6
#-从结果中提取网络节点和边数据
node = res[[2]][[2]]
head(node)## ID X1 X2 elements igraph.degree igraph.closeness
## 1 ASV_101 1.8922077 -16.17873 ASV_101 3 1
## 2 ASV_1050 10.4488008 11.60225 ASV_1050 1 1
## 3 ASV_113 -0.3111552 18.45045 ASV_113 0 0
## 4 ASV_1148 11.7082072 -13.80247 ASV_1148 0 0
## 5 ASV_1149 2.1328120 18.18372 ASV_1149 1 1
## 6 ASV_115 7.1706596 -16.75757 ASV_115 0 0
## igraph.betweenness igraph.cen.degree group Kingdom Phylum
## 1 0 3 B.T1.WT Bacteria Proteobacteria
## 2 0 1 B.T1.WT Bacteria Proteobacteria
## 3 0 0 B.T1.WT Bacteria Actinobacteria
## 4 0 0 B.T1.WT Bacteria Firmicutes
## 5 0 1 B.T1.WT Bacteria Firmicutes
## 6 0 0 B.T1.WT Bacteria Proteobacteria
## Class Order Family Genus
## 1 Unassigned Unassigned Unassigned Unassigned
## 2 Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas
## 3 Actinobacteria Actinomycetales Microbacteriaceae Agromyces
## 4 Bacilli Bacillales Planococcaceae Lysinibacillus
## 5 Bacilli Bacillales Paenibacillaceae_1 Unassigned
## 6 Betaproteobacteria Burkholderiales Oxalobacteraceae Unassigned
## Species Group time space links nodes
## 1 Unassigned WT T1 B 162 100
## 2 Sphingomonas_daechungensis WT T1 B 162 100
## 3 Unassigned WT T1 B 162 100
## 4 Unassigned WT T1 B 162 100
## 5 Unassigned WT T1 B 162 100
## 6 Unassigned WT T1 B 162 100
## label
## 1 B.T1.WT: (nodes: 100; links: 162)
## 2 B.T1.WT: (nodes: 100; links: 162)
## 3 B.T1.WT: (nodes: 100; links: 162)
## 4 B.T1.WT: (nodes: 100; links: 162)
## 5 B.T1.WT: (nodes: 100; links: 162)
## 6 B.T1.WT: (nodes: 100; links: 162)edge = res[[2]][[3]]
head(edge)## X2 Y2 OTU_2 OTU_1 weight X1 Y1 cor group Group
## 1 6.251566 6.418839 ASV_30 ASV_21 -1 5.453796 4.025147 - R.T3.OE OE
## 2 6.251566 6.418839 ASV_30 ASV_29 1 4.842470 7.579802 + R.T3.OE OE
## 3 6.251566 6.418839 ASV_30 ASV_55 1 3.675405 8.282087 + R.T3.OE OE
## 4 6.251566 6.418839 ASV_30 ASV_290 1 3.529926 5.996121 + R.T3.OE OE
## 5 4.842470 7.579802 ASV_29 ASV_21 -1 5.453796 4.025147 - R.T3.OE OE
## 6 4.842470 7.579802 ASV_29 ASV_30 1 6.251566 6.418839 + R.T3.OE OE
## time space links nodes label
## 1 T3 R 160 108 R.T3.OE: (nodes: 108; links: 160)
## 2 T3 R 160 108 R.T3.OE: (nodes: 108; links: 160)
## 3 T3 R 160 108 R.T3.OE: (nodes: 108; links: 160)
## 4 T3 R 160 108 R.T3.OE: (nodes: 108; links: 160)
## 5 T3 R 160 108 R.T3.OE: (nodes: 108; links: 160)
## 6 T3 R 160 108 R.T3.OE: (nodes: 108; links: 160)cb_palette <- c("#ed1299", "#09f9f5", "#246b93", "#cc8e12", "#d561dd", "#c93f00", "#ddd53e",
"#4aef7b", "#e86502", "#9ed84e", "#39ba30", "#6ad157", "#8249aa", "#99db27", "#e07233", "#ff523f",
"#ce2523", "#f7aa5d", "#cebb10", "#03827f", "#931635", "#373bbf", "#a1ce4c", "#ef3bb6", "#d66551",
"#1a918f", "#ff66fc", "#2927c4", "#7149af" ,"#57e559" ,"#8e3af4" ,"#f9a270" ,"#22547f", "#db5e92",
"#edd05e", "#6f25e8", "#0dbc21", "#280f7a", "#6373ed", "#5b910f" ,"#7b34c1" ,"#0cf29a" ,"#d80fc1",
"#dd27ce", "#07a301", "#167275", "#391c82", "#2baeb5","#925bea", "#63ff4f")
tem1 = node[,fill] %>% as.matrix() %>% as.vector() %>% unique()
tabf = data.frame(id = tem1,color =cb_palette[1:length(tem1)] )
head(tabf)## id color
## 1 Proteobacteria #ed1299
## 2 Actinobacteria #09f9f5
## 3 Firmicutes #246b93
## 4 Bacteroidetes #cc8e12
## 5 Unassigned #d561dd
## 6 Verrucomicrobia #c93f00colnames(tabf)[1] = fill
node1 = node %>% left_join(tabf,by = fill)
p2 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,
color = cor),
data = edge, size = 0.03,alpha = 0.5) +
geom_point(aes(X1, X2,
size = !!sym(size),
fill = color ),
pch = 21, data = node1,color = "gray40") +
facet_wrap(.~ label,scales="free_y",ncol = row.num ) +
# geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
# geom_text(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
scale_colour_manual(values = c("#6D98B5","#D48852")) +
scale_fill_manual(values = tabf$color,labels = tabf[,fill]) +
scale_size(range = c(0.8, maxnode)) +
scale_x_continuous(breaks = NULL) +
scale_y_continuous(breaks = NULL) +
theme(panel.background = element_blank(),
plot.title = element_text(hjust = 0.5)
) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()
) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p2
# ggsave("cs1.pdf",p2,width = 5*5,height = 4*3)
#-------另外一组网络绘制
data(ps)
res2 = Facet.network (
ps.st= ps,# phyloseq对象
g1 = "Group",# 分组1
# g2 = "space",# 分组2
# g3 = "time",# 分组3
ord.g1 = c("WT","KO","OE"),# 排序顺序
# ord.g2 = c("B","R") ,# 排序顺序
# ord.g3 = c("T1","T2","T3") ,# 排序顺序
order = "time", # 出图每行代表的变量
fill = "Phylum",
size = "igraph.degree",
layout_net = "model_maptree2",
r.threshold=0.8,
p.threshold=0.01,
method = "spearman",
select_layout = TRUE,
clu_method = "cluster_fast_greedy",
maxnode = 5
)## [1] "WT" "KO" "OE"p3 = res2[[1]]
p3
#-提取节点和边数据
node = res2[[2]][[2]]
head(node)## ID X1 X2 elements igraph.degree igraph.closeness
## 1 ASV_1 6.394069 -0.02808431 ASV_1 1 1
## 2 ASV_10 3.652617 -21.02187428 ASV_10 0 0
## 3 ASV_100 7.078794 10.35742308 ASV_100 1 1
## 4 ASV_101 -3.792132 16.03963922 ASV_101 1 1
## 5 ASV_102 8.292658 -9.25473746 ASV_102 1 1
## 6 ASV_103 4.025401 21.04072916 ASV_103 0 0
## igraph.betweenness igraph.cen.degree group Kingdom Phylum
## 1 0 1 ..WT Bacteria Actinobacteria
## 2 0 0 ..WT Bacteria Proteobacteria
## 3 0 1 ..WT Bacteria Proteobacteria
## 4 0 1 ..WT Bacteria Proteobacteria
## 5 0 1 ..WT Bacteria Proteobacteria
## 6 0 0 ..WT Bacteria Proteobacteria
## Class Order Family Genus
## 1 Actinobacteria Actinomycetales Thermomonosporaceae Unassigned
## 2 Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium
## 3 Betaproteobacteria Burkholderiales Comamonadaceae Hydrogenophaga
## 4 Unassigned Unassigned Unassigned Unassigned
## 5 Gammaproteobacteria Unassigned Unassigned Unassigned
## 6 Alphaproteobacteria Rhizobiales Phyllobacteriaceae Phyllobacterium
## Species Group time space links nodes
## 1 Unassigned WT 138 105
## 2 Unassigned WT 138 105
## 3 Hydrogenophaga_intermedia WT 138 105
## 4 Unassigned WT 138 105
## 5 Unassigned WT 138 105
## 6 Phyllobacterium_bourgognense WT 138 105
## label
## 1 ..WT: (nodes: 105; links: 138)
## 2 ..WT: (nodes: 105; links: 138)
## 3 ..WT: (nodes: 105; links: 138)
## 4 ..WT: (nodes: 105; links: 138)
## 5 ..WT: (nodes: 105; links: 138)
## 6 ..WT: (nodes: 105; links: 138)edge = res2[[2]][[3]]
head(edge)## X2 Y2 OTU_2 OTU_1 weight X1 Y1 cor group
## 1 1.9872552 1.5012669 ASV_22 ASV_32 1 -3.0932354 -2.3169033 + ..OE
## 2 1.9872552 1.5012669 ASV_22 ASV_77 1 -0.4037329 -0.8876355 + ..OE
## 3 1.9872552 1.5012669 ASV_22 ASV_95 -1 -0.9859822 2.3370501 - ..OE
## 4 1.9872552 1.5012669 ASV_22 ASV_214 1 2.3530761 -1.0455685 + ..OE
## 5 1.9872552 1.5012669 ASV_22 ASV_156 1 -3.4791850 0.5142463 + ..OE
## 6 -0.4037329 -0.8876355 ASV_77 ASV_32 1 -3.0932354 -2.3169033 + ..OE
## Group time space links nodes label
## 1 OE 148 96 ..OE: (nodes: 96; links: 148)
## 2 OE 148 96 ..OE: (nodes: 96; links: 148)
## 3 OE 148 96 ..OE: (nodes: 96; links: 148)
## 4 OE 148 96 ..OE: (nodes: 96; links: 148)
## 5 OE 148 96 ..OE: (nodes: 96; links: 148)
## 6 OE 148 96 ..OE: (nodes: 96; links: 148)node$group %>% unique()## [1] "..WT" "..KO" "..OE"#----第二组网络用第一组颜色填充#---------
tem1 = node[,fill] %>% as.matrix() %>% as.vector() %>% unique()
tabf2 = data.frame(id = tem1 )
colnames(tabf2)[1] = fill
head(tabf2)## Phylum
## 1 Actinobacteria
## 2 Proteobacteria
## 3 Bacteroidetes
## 4 Unassigned
## 5 Chloroflexi
## 6 Firmicutes# tabf2[1,1] = "wentao"
tabf3 = tabf2 %>% left_join(tabf,by = fill)
num = tabf3$color[is.na(tabf3$color)] %>% length()
if (num == 0) {
tabf.f = tabf3
}else{
tem = tabf3 %>% filter(is.na(color))
tem$color = setdiff(cb_palette,tabf$color)[1:length(tem$color)]
tabf.f = rbind(tabf3,tem)
}
colnames(tabf.f)[1] = fill
node1 = node %>% left_join(tabf.f,by = fill)
head(node1)## ID X1 X2 elements igraph.degree igraph.closeness
## 1 ASV_1 6.394069 -0.02808431 ASV_1 1 1
## 2 ASV_10 3.652617 -21.02187428 ASV_10 0 0
## 3 ASV_100 7.078794 10.35742308 ASV_100 1 1
## 4 ASV_101 -3.792132 16.03963922 ASV_101 1 1
## 5 ASV_102 8.292658 -9.25473746 ASV_102 1 1
## 6 ASV_103 4.025401 21.04072916 ASV_103 0 0
## igraph.betweenness igraph.cen.degree group Kingdom Phylum
## 1 0 1 ..WT Bacteria Actinobacteria
## 2 0 0 ..WT Bacteria Proteobacteria
## 3 0 1 ..WT Bacteria Proteobacteria
## 4 0 1 ..WT Bacteria Proteobacteria
## 5 0 1 ..WT Bacteria Proteobacteria
## 6 0 0 ..WT Bacteria Proteobacteria
## Class Order Family Genus
## 1 Actinobacteria Actinomycetales Thermomonosporaceae Unassigned
## 2 Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium
## 3 Betaproteobacteria Burkholderiales Comamonadaceae Hydrogenophaga
## 4 Unassigned Unassigned Unassigned Unassigned
## 5 Gammaproteobacteria Unassigned Unassigned Unassigned
## 6 Alphaproteobacteria Rhizobiales Phyllobacteriaceae Phyllobacterium
## Species Group time space links nodes
## 1 Unassigned WT 138 105
## 2 Unassigned WT 138 105
## 3 Hydrogenophaga_intermedia WT 138 105
## 4 Unassigned WT 138 105
## 5 Unassigned WT 138 105
## 6 Phyllobacterium_bourgognense WT 138 105
## label color
## 1 ..WT: (nodes: 105; links: 138) #09f9f5
## 2 ..WT: (nodes: 105; links: 138) #ed1299
## 3 ..WT: (nodes: 105; links: 138) #ed1299
## 4 ..WT: (nodes: 105; links: 138) #ed1299
## 5 ..WT: (nodes: 105; links: 138) #ed1299
## 6 ..WT: (nodes: 105; links: 138) #ed1299p4 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,
color = cor),
data = edge, size = 0.03,alpha = 0.5) +
geom_point(aes(X1, X2,
fill =color,
size = !!sym(size)),
pch = 21, data = node1,color = "gray40") +
facet_wrap(.~ label,scales="free_y",ncol = row.num ) +
# geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
# geom_text(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
scale_colour_manual(values = c("#6D98B5","#D48852")) +
scale_fill_manual(values = tabf.f$color,labels = tabf.f[,fill]) +
scale_size(range = c(0.8, maxnode)) +
scale_x_continuous(breaks = NULL) +
scale_y_continuous(breaks = NULL) +
theme(panel.background = element_blank(),
plot.title = element_text(hjust = 0.5)
) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()
) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
# guides(fill=FALSE)
p4
# ggsave("cs2.pdf",p4,width = 5*3,height = 4*1)library(patchwork)
p2|p4